The present invention relates to the numerical simulation of fluid flow and, more particularly, to a numerical method of simulating barotropic fluid flow past a body by minimizing a potential functional.
Barotropic fluid flow is described by the Euler equations:                                                         ∂                              v                →                                                    ∂              t                                +                                    v              →                        ·                                          ∇                →                            ⁢                              v                →                                              +                                    ∇              →                        ⁢                          (                              h                +                Φ                            )                                      =        0                            (        1        )            
and of the continuity equation                                                         ∂              ρ                                      ∂              t                                +                                    ∇              →                        ·                          (                              ρ                ⁢                                  xe2x80x83                                ⁢                                  v                  →                                            )                                      =        0                            (        2        )            
where {right arrow over (xcexd)}(xk,t) is the velocity vector field, "PHgr"(xk,t) is a potential relating to an external force, xcfx81(xk,t) is the scalar density field, and h, the specific enthalpy, is a function of xcfx81 through the equation of state of the fluid. {right arrow over (xcexd)}, "PHgr" and xcfx81 are functions of position xk=(x,y,z) in space and of time t. In the special case of a stationary flow with small viscosity, the time derivatives vanish:
{right arrow over (xcexd)}xc2x7{right arrow over (∇)}{right arrow over (xcexd)}+{right arrow over (∇)}(h+"PHgr")=0; {right arrow over (∇)}xc2x7(xcfx81{right arrow over (xcexd)})=0xe2x80x83xe2x80x83(3)
Only in trivial cases can these equations be solved analytically. In cases of practical interest, for example, in aerodynamics and hydrodynamics, these equations must be solved numerically. The most common way to solve these equations is by integrating them. The possible numerical instabilities associated with numerical integration are well known and need not be detailed here. A numerical solution based on a variational principle would be inherently numerically stable. Ideally, such a numerical solution would involve finding the minima of a potential functional with respect to the associated variables. In the case of fluid flow, the variables are four scalar fields: xcfx81 and the three Cartesian components of {right arrow over (xcexd)}. Heretofore, no variational formulation of barotropic fluid flow has represented the solutions of equations (1) and (2) or of equations (3) as the values of {right arrow over (xcexd)} and xcfx81, or of any other set of only four scalar fields, that minimize a potential functional. The smallest number of scalar fields obtained heretofore was seven, by R. L. Seliger and G. B. Witham (Proc. Roy. Soc. London, Vol. A305, p. 1, 1968). Because equations (1)-(3) depend on only four scalar fields, it should be possible to obtain a variational formulation of barotropic fluid flow in terms of only four scalar fields.
There is thus a widely recognized need for, and it would be highly advantageous to have, a method of finding numerical solutions of the equations of fluid flow by minimizing a potential functional of four scalar fields.
According to the present invention there is provided a numerical method of simulating fluid flow past a body having a boundary, including the steps of: (a) formulating the fluid flow in terms of a potential functional of at most four fundamental scalar fields and at most two scalar variables; and (b) extremizing the potential functional with respect to the at most four fundamental scalar fields and with respect to the at most two scalar variables.
According to the present invention there is provided a system for numerical simulation of fluid flow, including: (a) a software module including a plurality of instructions for computing, from discrete representations of at most four fundamental scalar fields, a potential functional representative of the fluid flow, and for varying the discrete representations to extremize the potential functional; (b) a processor for executing the instructions; and (c) a memory for storing the discrete representations.
According to the present invention there is provided a numerical method of simulating fluid flow, including the steps of: (a) formulating the fluid flow in terms of a potential functional of at most four fundamental scalar fields and at most two scalar variables; and (b) extremizing the potential functional with respect to the at most four fundamental scalar fields and with respect to the at most two scalar variables.
According to the present invention there is provided a numerical method of simulating fluid flow past a body having a boundary, including the steps of: (a) formulating the fluid flow in terms of an action which is the sum of: (i) an action integral, from an initial time t0 to a final time t1, of a Lagrangian functional L of at most four fundamental scalar fields, (ii) a spatial integral of a function of the at most four fundamental scalar fields at the final time t1, and (iii) a negative of a spatial integral of the function of the at most four fundamental scalar fields at the initial time t0; and (b) extremizing the action with respect to the at most four fundamental scalar fields.
According to the present invention there is provided a numerical method of simulating fluid flow including the steps of: (a) formulating the fluid flow in terms of an action which is the sum of: (i) an action integral, from an initial time t0 to a final time t1, of a Lagrangian functional L of at most four fundamental scalar fields, (ii) a spatial integral of a function of the at most four fundamental scalar fields at the final time t1, and (iii) a negative of a spatial integral of the function of the at most four fundamental scalar fields at the initial time t0; and (b) extremizing the action with respect to the at most four fundamental scalar fields.
According to the present invention there is provided a system for numerical simulation of fluid flow, including: (a) a software module including a plurality of instructions for computing, from discrete representations of at most four fundamental scalar fields, an action representative of the fluid flow, and for varying the discrete representations to extremize the action; (b) a processor for executing the instructions; and (c) a memory for storing the discrete representations.
It is shown in the Theory section below that stationary barotropic fluid flow can be formulated as the extremization of a potential functional of at most four scalar fields: the Clebsch scalar fields xcex1, xcex2 and xcexd and the density field xcfx81; and also of at most two scalar variables xcex21 and xcexd1. Specifically, the functional is minimized with respect to xcex1, xcex2 and xcexd and maximized with respect to xcex21 and xcexd1. The scalar fields are functions of the position vector xk. The velocity vector field {right arrow over (xcexd)} is related to the Clebsch scalar fields through {right arrow over (xcexd)}=xcex1{right arrow over (∇)}xcex2+{right arrow over (∇)}xcexd.
The potential functionals are as follows:
If the component of {right arrow over (xcexd)} normal to the boundary of the region of fluid flow is zero, or if the density xcfx81 on the boundary is zero, the potential functional is                                           ∫            V                    ⁢                                    [                                                                    1                    2                                    ⁢                                                            (                                                                        α                          ⁢                                                      xe2x80x83                                                    ⁢                                                                                    ∇                              →                                                        ⁢                            β                                                                          +                                                                              ∇                            →                                                    ⁢                          v                                                                    )                                        2                                                  +                                  ϵ                  ⁢                                      xe2x80x83                                    ⁢                                      (                    ρ                    )                                                  +                Φ                            ]                        ⁢            ρ            ⁢                          xe2x80x83                        ⁢                          ⅆ                                                                   3                                ⁢                x                                                    -                              β            1                    ⁡                      (                                                            ∫                  V                                ⁢                                  α                  ⁢                                      xe2x80x83                                    ⁢                                      ρ                    0                                    ⁢                                      xe2x80x83                                    ⁢                                      ⅆ                                                                                           3                                            ⁢                      x                                                                                  -                                                α                  _                                ⁢                                  xe2x80x83                                ⁢                                  M                  0                                                      )                          -                              v            1                    ⁡                      (                                                            ∫                  V                                ⁢                                  ρ                  ⁢                                      xe2x80x83                                    ⁢                                      ⅆ                                                                                           3                                            ⁢                      x                                                                                  -                              M                0                                      )                                              (        4        )            
(equations 56 and 62 of the Theory section below) for compressible flows and                                           ∫            V                    ⁢                                    [                                                                    1                    2                                    ⁢                                                            (                                                                        α                          ⁢                                                      xe2x80x83                                                    ⁢                                                                                    ∇                              →                                                        ⁢                            β                                                                          +                                                                              ∇                            →                                                    ⁢                          v                                                                    )                                        2                                                  +                Φ                            ]                        ⁢                          ρ              0                        ⁢                          xe2x80x83                        ⁢                          ⅆ                                                                   3                                ⁢                x                                                    -                              β            1                    ⁡                      (                                                            ∫                  V                                ⁢                                  α                  ⁢                                      xe2x80x83                                    ⁢                                      ρ                    0                                    ⁢                                      xe2x80x83                                    ⁢                                      ⅆ                                                                                           3                                            ⁢                      x                                                                                  -                                                α                  _                                ⁢                                  xe2x80x83                                ⁢                                  M                  0                                                      )                                              (        5        )            
(equations 65 and 66 of the Theory section below) for incompressible flows. xcex5(xcfx81) is the specific internal energy of the fluid, in the case of a compressible flow, and is related to the specific enthalpy by.                     h        =                              ∂                          (                              ρ                ⁢                                  xe2x80x83                                ⁢                ϵ                            )                                            ∂            ρ                                              (        6        )            
{overscore (xcex1)} and M0 are constants. The integrals are over the volume V occupied by the fluid. xcfx810 is the constant density, in the case of an incompressible flow. The functional is minimized with respect to xcex1, xcex2 and xcexd, and maximized with respect to xcex21. In the case of a compressible flow, the functional also is minimized with respect to xcfx81 and maximized with respect to xcexd1.
In the case of potential flows, for which {right arrow over (∇)}xc3x97{right arrow over (xcexd)}=0, the potential functional is                                           ∫            V                    ⁢                                    [                                                                    1                    2                                    ⁢                                                            (                                                                        ∇                          →                                                ⁢                                                  v                          ^                                                                    )                                        2                                                  +                                  ϵ                  ⁢                                      xe2x80x83                                    ⁢                                      (                    ρ                    )                                                  +                Φ                            ]                        ⁢            ρ            ⁢                          xe2x80x83                        ⁢                          ⅆ                                                                   3                                ⁢                x                                                    -                              v            1                    ⁡                      (                                                            ∫                  V                                ⁢                                  ρ                  ⁢                                      xe2x80x83                                    ⁢                                      ⅆ                                                                                           3                                            ⁢                      x                                                                                  -                              M                0                                      )                                              (        7        )            
(equations 102 and 103 of the Theory section below) for compressible flows and                               ∫          V                ⁢                              [                                                            1                  2                                ⁢                                                      (                                                                  ∇                        →                                            ⁢                      v                                        )                                    2                                            +              Φ                        ]                    ⁢                      ρ            0                    ⁢                      xe2x80x83                    ⁢                      ⅆ                                                           3                            ⁢              x                                                          (        8        )            
(equation 106 of the Theory section below) for incompressible flows. xcex5(xcfx81) is the specific internal energy of the fluid, in the case of a compressible flow. The integrals are over the volume occupied by the fluid. M0 is a constant. xcfx810 is the constant density, in the case of an incompressible flow. The functional is minimized with respect to xcexd. In the case of a compressible flow, the functional also is minimized with respect to xcfx81 and maximized with respect to xcexd1.
It also is shown in the Theory section below (see equation 113 therein) that non-stationary barotropic fluid flow can be formulated as the maximization of an action which is the sum of three integrals: an action integral of a Lagrangian functional L of three Clebsch scalar fields xcex1, xcex2 and xcexd and the density field xcfx81 from an initial time t0 to a final time t1.                               ∫                      t            0                                t            1                          ⁢                  L          ⁢                      xe2x80x83                    ⁢                      ⅆ            t                                              (        9        )            
an initial integral over the four scalar fields at time t0:                     -                              ∫            V                    ⁢                                    (                                                                    α                    ⁡                                          (                                                                        x                          k                                                ,                                                  t                          0                                                                    )                                                        ⁢                                      β                    ⁡                                          (                                                                        x                          k                                                ,                                                  t                          0                                                                    )                                                                      +                                  v                  ⁡                                      (                                                                  x                        k                                            ,                                              t                        0                                                              )                                                              )                        ⁢                          ρ              ⁡                              (                                                      x                    k                                    ,                                      t                    0                                                  )                                      ⁢                          xe2x80x83                        ⁢                          ⅆ                                                                   3                                ⁢                x                                                                        (        10        )            
and a similar final integral over the four scalar fields at time t1:                               ∫          V                ⁢                              (                                                            α                  ⁡                                      (                                                                  x                        k                                            ,                                              t                        1                                                              )                                                  ⁢                                  β                  ⁡                                      (                                                                  x                        k                                            ,                                              t                        1                                                              )                                                              +                              v                ⁡                                  (                                                            x                      k                                        ,                                          t                      1                                                        )                                                      )                    ⁢                      ρ            ⁡                          (                                                x                  k                                ,                                  t                  1                                            )                                ⁢                      xe2x80x83                    ⁢                      ⅆ                                                           3                            ⁢              x                                                          (        11        )            
The spatial integrals are over the volume V occupied by the fluid. Note that the integrand of the initial integral is the negative of the integrand of the final integral.
If the component of {right arrow over (xcexd)} normal to the boundary of the region of fluid flow is zero, or if the density xcfx81 on the boundary is zero, the Lagrangian functional is:                               ∫          V                ⁢                              [                                          -                                  (                                                            α                      ⁢                                              xe2x80x83                                            ⁢                                                                        ∂                          β                                                                          ∂                          t                                                                                      +                                                                  ∂                        v                                                                    ∂                        t                                                                              )                                            -                                                1                  2                                ⁢                                  xe2x80x83                                ⁢                                                      (                                                                  α                        ⁢                                                                              ∇                            →                                                    ⁢                          β                                                                    +                                                                        ∇                          →                                                ⁢                        v                                                              )                                    2                                            -                              ϵ                ⁢                                  xe2x80x83                                ⁢                                  (                  ρ                  )                                            -              Φ                        ]                    ⁢          ρ          ⁢                      xe2x80x83                    ⁢                      ⅆ                                                           3                            ⁢              x                                                          (        12        )            
(equation 31 of the Theory section below) for compressible flows and                               ∫          V                ⁢                              [                                          -                                  (                                                            α                      ⁢                                              xe2x80x83                                            ⁢                                                                        ∂                          β                                                                          ∂                          t                                                                                      +                                                                  ∂                        v                                                                    ∂                        t                                                                              )                                            -                                                1                  2                                ⁢                                  xe2x80x83                                ⁢                                                      (                                                                  α                        ⁢                                                                              ∇                            →                                                    ⁢                          β                                                                    +                                                                        ∇                          →                                                ⁢                        v                                                              )                                    2                                            -              Φ                        ]                    ⁢                      ρ            0                    ⁢                      xe2x80x83                    ⁢                      ⅆ                                                           3                            ⁢              x                                                          (        12        )            
(equation 35 of the Theory section below) for incompressible flows. As in the case of stationary flows, xcex5(xcfx81) is the specific internal energy of the fluid, in the case of a compressible flow; the integrals are over the volume V occupied by the fluid; and xcfx810 is the constant density, in the case of an incompressible flow. The maximization is with respect to xcex1, xcex2 and xcexd, and also with respect to xcfx81 in the case of a compressible flow.
The Clebsch scalar fields xcex1, xcex2 and xcexd and the density xcfx81 are referred to herein as xe2x80x9cfundamentalxe2x80x9d scalar fields because only these fields are varied explicitly to minimize the potential functionals. The other scalar fields that appear in the potential functionals are functions of these fundamental scalar fields. For example, xcex5 is a function of xcfx81, via equation (6).
In addition to being inherently stable, the present invention, as applied to stationary flows, is faster than the prior art methods and is easier to use for complex geometries.
A Note on Notation
In the Theory section below, the symbols xcex1, xcex2 and xcexd are used to represent Clebsch variables that are functions of both space xk and time t. In the case of stationary flows, the space and time dependences of the Clebsch variables can be separated as follows:
xcex1={circumflex over (xcex1)}(xk)xe2x80x83xe2x80x83(14)
xcex2={circumflex over (xcex2)}(xk)xe2x88x92xcex21txe2x80x83xe2x80x83(15)
xcexd={circumflex over (xcexd)}(xk)xe2x88x92xcexd1txe2x80x83xe2x80x83(16)
In potential functionals for stationary flows as derived in the Theory section below, the Clebsch scalar fields are these functions {circumflex over (xcex1)}, {circumflex over (xcex2)} and {circumflex over (xcexd)} of xk only; and in the course of the derivations, the dependence on time drops out, leaving xcex21 and xcexd1, as scalar variables to be varied, along with the values of the scalar fields {circumflex over (xcex1)}, {circumflex over (xcex2)} and {circumflex over (xcexd)}, to extremize the potential functionals. For simplicity elsewhere herein, when the potential functionals of stationary flows are discussed, the Clebsch scalar fields that enter into these potential functionals are denoted by xcex1, xcex2 and xcexd, without carets.
The symbol that is used herein for velocity, a lowercase italic xe2x80x9cVxe2x80x9d (xcexd), and the symbol that is used herein for the third Clebsch scalar field, a lowercase Greek xe2x80x9cnuxe2x80x9d (xcexd), are typographically similar. To avoid confusion, the reader should bear in mind that the symbol for velocity always appears with an arrow above it ({right arrow over (xcexd)}) when referring to the vectorial velocity field, or with the subscripts x or y when referring to Cartesian components of the velocity field; and that the symbol for the third Clebsch scalar field appears either unadorned (xcexd) when referring to the third Clebsch scalar field itself, or with the subscript 1 when referring to a constant that is related to the third Clebsch scalar field.